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1. 


INTRODUCTION 


Current  Air  Force  interest  is  high  on  ramjet  technology  due 
to  the  advent  of  the  National  Aerospace  Plane  project  (NASP) . 

The  aerodynamic  community  is  experiencing  an  exponential  growth 
in  Computational  Fluid  Dynamics  ( CFD)  as  witnessed  by  the  large 
number  of  publications  now  appearing  in  this  field.  It  is, 
therefore,  appropriate  that  CFD  methodology  be  brought  to  bear  on 
the  ramjet  technology.  By  numerical  modeling  the  flow  process 
and  conducting  validating  experiments  the  ramjet  field  may  be 
advanced.  This  report  documents  work  on  the  Aero  Propulsion 
Laboratory  Contract  (No.  F33615-86-C-2615)  on  the  computation  of 
Ramjet  Internal  Flowfields. 

2.  R&D  STATUS  REPORT 

The  CFD  effort  on  Ramjet  Internal  Flowfields  may  be  divided 
into  two  categories,  i.e.  subsonic  combustion  and  supersonic 
combustion.  The  first  category  (subsonic)  employs  the  time- 
dependent  Navier-Stokes  code  and  adds  the  species  equations  for 
hydrocarbon  fuel  and  air  combustion.  A  forward  reaction  rate 
model  is  used  for  the  complete  chemical  reaction.  The  second 
category  (supersonic)  utilizes  the  parabolized  Navier-Stokes 
program  with  additional  species  equations  for  hydrogen/air 
combustion.  The  "extent  of  reaction"  method  is  used  to  account 
for  equilibrium  chemistry.  These  two  categories  now  will  be 
discussed. 

2 . 1  Subsonic  Combustion 

2.1.1  Governing  Equations 

The  governing  equations  for  combustion  driven 
flows  include  the  conservation  of  species,  momentum  and  energy, 
plus  the  state  equation.  These  are  listed  below  in  conservation 
form: 

3pY. 

IT  +  V*^Y iY  -  PD.jVY.)  =  m. 

dpV 

+  V. (pvv  -  r )  =  o 


l 


dge 

at 


V. (pVe  -  L'Y  -  g) 


0 


P  -  £  Pi  -  PRT  t  ^  ;  h  =  I  Yihi 


1 


where 


1  =  -  (P  +  r  mV  -V)  I  +  M(VV  +  VV ) 


e  =  h  +  ?-  -  B 

2  p 


g  =  kVT  +  £  MDiihiVYi 


ID 

MC 


=  S 


pDij  °c.  ' 


k  =  Pr 


mi  = 


hi  =  J  CpidT  +  hi 

\f  ,  1/  . 

-k.  e  Ei/RT  (p  x) (p  Dj  Arrhenius  Law 

*■  J 


mj  =  Sjin1  ,  Sj  =  Stoichiometric  Coefficient 


2.1.2  Transport  Properties  of  a  Mixture 

To  solve  the  conservation  equations  for  a 
mixture  of  gases  the  transport  coefficients  (p,  k,  D^)  are 
required.  A  method  for  determining  these  coefficients  has  been 
reported  by  Wilke1,  however,  for  the  subsonic  ramjet  in  which  75 
percent  of  the  mixture  is  nitrogen,  the  following,  simpler, 
procedure  will  be  used. 

Sutherland's  law  for  viscosity: 

*  =  ^(rVi) 


2 


-8  lb  sec 

where  p.  =  2.27  x  10 

r  ftc 


air  values 

i 


S  =  198. 6°R 


M _ 


S 


c 


i 


2 

For  T  =  1 , 800 *  R,  Kanury  presents  the  following  values  for  air. 

Pr  =  0.70 

sc  =  1.00 

S 

(Note  this  corresponds  to  a  Lewis  number,  Le  =  ^  =  1.4.)  These 

r 

values  appear  appropriate  for  the  early  calculations.  The  values 
for  the  molecular  transport  phenomena  are  masked  by  the  unsteady 
convection  transport  being  simulated  by  the  time-dependent 
oscillations.  Therefore,  the  computation  is  relatively 
insensitive  to  the  values  of  these  molecular  transport 
coefficients . 


2.1.3  Thermodynamic  Properties 

In  the  CFD  codes  involving  different  gas 
species  at  elevated  temperatures  there  is  a  need  to  evaluate  the 
thermodynamic  properties  as  a  function  of  temperature.  The 
procedure  should  be  as  simple  as  possible  and  not  require 
iteration. 


The  conservation  equations  provide  the 
following  information  at  each  point  in  time  and  space. 

p,  p u,  pv,  p w,  pe,  p Yi 


3 


From  this  information  it  is  necessary  to  determine  p  and  T. 
following  equations  provide  this  information: 

P  =  z  Pi  Where  pi  =  —  pRT 


where  R  =  1.987  cal/(gm  mole) ‘K  or  49,723 


2 

sec  °R  mole 


=  molecular  weight 


Hence  p  =  p RT  where  R  =  [  ^  R. 


The  temperature  must  be  obtained  from 


enthalpy. 


“  I  Yi  hi 
rT 

=  C  dT  +  h? 

Jo  Pi  1 


(based  on  mass) 


C  =  specific  heat  at  constant  pressure  (based  on  mole) 


Therefore, 


h  -  C  T  +  l  Y.h“ 


where 


Yi  rT 

c  =  Y  C  dT 

P  L  T  J0  p. 


Hence,  temperature  may  be  deduced  from  the  known  value  of 
internal  energy  (e) . 


e  =  h.B  +  Yl  = 
e  p  2 


CpT  +  l  Y-hV  -  RT  +  2 
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e  - 


or 


T  = 


v_ 

2 


-  I  Y.h? 
L  11 


-  R 


The  complexity  in  this  equation  is  due  to  the  fact  that  is  a 
function  of  T,  requiring  an  iteration  to  resolve.  Several 
methods  of  solution  are  possible,  each  with  different 
computational  efficiency.  The  simplest  method  is  to  evaluate 
at  time  (t  -  At)  to  find  T't).  The  alternate  approach  is  to 
solve  a  low  order  polynomial  for  the  positive  real  root.  To 
obtain  C  (T)  a  quadratic  polynomial  for  the 

C  for  each  species  was  derived.  Using  the  JANNAF  Tables  of 

pi 

Thermochemical  Data,  a  best  fit  for  the  five  species  involved  in 
the  combustion  of  JP-4  and  air  was  accomplished  over  the 
temperature  range  from  900°  to  4,500°R. 

C  =a.  +b.T+c.T2 

p.  l  l  l 


These  values  are  tabulated  in  Table  1.  The  global  value  for  C 

P 

may  be  obtained  by  the  following  integration: 


Yi  fT  T 

=  l  —  C  dT  =  l  Y .  a .  +  \  l 

J0pi 


T 

Yibi  +  r 


Yici 


Each  species  will  require  seven  quantities  to 
determine  all  property  values.  For  hydrocarbon  combustion  over  a 
temperature  range  from  900°  to  4,500°R  the  values  in  Table  1  are 
appropriate . 
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TABLE  1 

SUMMARY  OF  JANNAF  TABLES  OF  THERMOCHEMICAL  DATA 


Species 

M. 

l 

CP  -  ai  +  biT 

ft2/sec2 

2 

+  CiT 

•R 

2  ,  2 
ft  /sec 

*hi 

m  .  =  S  . ; 
D  D 

S  . 

l 

bi 

Ci 

CH  2 

14 

10,740 

. 

6.396 

-7 . 4  6x 10~4 

-0. 160X108 

1.00 

°2 

32 

5,156 

0.927 

-1.09xl0~4 

0 

+3.429 

C02 

44 

3,623 

2.913 

-4 .48xl0-4 

-0. 964X108 

-3.143 

H2° 

18 

10,342 

1.856 

0 

-1.450X108 

-1.286 

N2 

28 

5,997 

0.517 

0 

0 

0 

R  =  49,732  ft2/sec2  °R  mole,  Pr  =  0.70, 
and  Rt  -  S- 


♦This  corresponds  to 
18.9  kBtu/lbm  fuel. 


2.1.4  Reaction  Rate  Model 

As  stated  previously,  an  Arrhenius  Law 
reaction  model  was  used  for  the  combustion  of  a  hydrocarbon. 
Locally  the  following  chemical  equation  prevails. 

CH2  +  1.5  02  -  C02  +  H20 

The  following  reaction  rate  model  is  used  to 
describe  the  reduction  of  the  hydrocarbon. 


where 


—  =  2 6 , 000° R  activation  temperature 


K 


k,  =  iO18/!41^)1-5  sec 


ft' 


The  remaining  species  production  will  be 
determined  from  the  stoichiometric  coefficients  (S^)  and  m^. 

m.  = 

The  values  of  are  listed  in  Table  1. 

2.1.5  Initial  Values 

Two  time-dependent  computations  for 
combustion  in  the  model  ramjet  were  accomplished.  The  initial 
conditions  for  the  cases  were:  TQ  =  1,000°R,  PQ  =  33  psia  at 
fuel-air  ratios  of  .020  and  .032.  The  mixture  entering  the 
combustor  is  vitiated  air  which  has  been  heated  to  1,000'R  by 
burning  hydrocarbon  fuel.  The  nitrogen/oxygen  ratio  is  returned 
to  3.76  (by  volume)  by  adding  02 . 

A  summary  of  the  mass  fractions,  for  the  five 
species  of  interest,  entering  the  combustor  region  for  the  two 
cases  is  shown  in  Table  2. 


The  upstream  boundary  conditions  are: 

PQ  =  33  psia  =  4,752  psf 
TQ  =  1,000'R 


v  =  0 

^  =  0 
dx 
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TABLE  2 

SUMMARY  OF  MASS  FRACTIONS 
ENTERING  THE  COMBUSTOR  REGION 


Species 

<t>  f  =  0.020 

l 

=  0.032 

l 

CH2 

0.020 

0.032 

°2 

0.222 

0.219 

C02 

0.018 

0.018 

H2° 

0.007 

0.007 

»2 

0.733 

0.724 

£  Yi  =  1.000  l  Yi  =  1.000 
2.1.6  Equilibrium  Chemistry 

Before  devising  a  scheme  to  represent  a 
combustion  model  it  is  necessary  to  examine  the  equilibrium  state 
at  the  anticipated  final  condition.  In  this  manner  the  appropriate 
species  may  be  identified. 


The  classical  procedure  involves  the  use  of  mole 

fractions  and  was  used  in  this  analysis,  even  though  the  CFD  code 

employs  mass  fractions.  The  mole  fractions  of  the  equilibrium 

composition  of  the  products  of  combustion  of  (CH  )  are  determined 

k 

using  the  equations  involving  equilibrium  constants  with  the 
assumption  that  the  component  gases  obey  the  equation  of  state  for 
an  ideal  gas. 


P 


i 


RT 


For  the  combustion  of  a  hydrocarbon  fuel  at  various  equivalence 
ratios  the  general  chemical  reaction  involving  ten  species  for 
products  may  be  written  as  follows: 

k  <CH2>  +  ^<°2  +  3‘76  N2>  * 


n 


N2<»2> 


+  n 


02<°2> 


+  nC02(C°2>  +  nH20(H2°»  +  nNO(NO> 


8 


9 


where  KD  .  =  f(T)  from  JANNAF  Tables.  Two  additional  variables  are 
introduced  in  these  equations,  i.e.,  p  and  n.  p  =  pressure  =  £  p. 
and  must  be  specified,  n  =  £  n^  =  total  number  of  moles  of  the 
combustion  products  and  becomes  an  auxiliary  equation. 

The  problem  has  now  been  formulated.  Given  0,  p, 
and  T  there  exists  11  algebraic  equations  involving  10  molar  values 
(n^)  plus  the  total  number  of  moles  (n) .  Conceptually ,  this  system 
of  equations  could  be  reduced  to  a  single  very  high  degree 
polynomial.  Physical  reasoning  requires  that  only  one  positive 
real  root  exist.  In  practice  the  system  is  reduced  to  fewer 
equations  by  substitution. 

The  following  iteration  sequence  was  used  to 
resolve  the  current  problem: 


n  =  y  n. 
u  1 


An  example  problem  has  been  solved  for  0  =  0.6  and  p  =  2 
atmospheres  at  a  temperature  range  of  900°  to  4,500°R.  This  state 
approximates  the  ramjet  combustor  condition.  A  summary  of  the 
result  is  shown  in  Tables  3  and  4. 

An  analysis  of  equilibrium  chemistry  for 
hydrocarbon  combustion  in  a  ramjet  has  been  accomplished.  It  was 
demonstrated  that  all  products  of  combustion  beyond  the  basic 
species  of  N2,  02,  H20,  and  C02  occur  at  traces  less  than  one 
percent  for  temperatures  below  4,140°F.  Based  upon  this  analysis 
it  appears  satisfactory  to  use  a  single-step  global  reaction  of  the 
following  form*: 

k<CH2>  +Lf  {°2  +  3'76  N2>  * 

fv 

(1^)3.76  N2  +  1.5(0  -  1)02  +  C02  +  H20 


* 

Note: 


0  <  1  for  a  lean  mixture. 
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TABLE  3 

HYDROCARBON-AIR  COMBUSTION 


.6  Equilibrium  States  p  =  2  atmospheres 


,500'  R 


1,083  38.9 


Jn 

yn 

4.903 

237 

yiT 

yiT 

0.0200 

0.0593 

0.558 

0.871 

0.00115yn 

0.0177yn 

4.7  x  10~4yn 

0.0102yn 

9.370 

9.313 

0.992 

0.931 

0.967 

0.923 

0.997 

0.914 

0.0602 

0.174 

0.0148 

0.101 

0.00330 

0.0864 

0.00163 

0.0347 

7.25  x  10-4 

0.0145 

1.09  x  10~4 

0.0075 

12.407 

12.499 
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TABLE  4 

MOLE  FRACTIONS 


<f>  =  0.6 

_ 2i _ 

P  =  2  atmospheres 

Species 

T  =  2 , 700  °  R 

T  =  3 , 600  °  R 

T  =  4,500-R 

M. 

1 

N2 

0.758 

0.755 

0.745 

28 

H2° 

0.0806 

0.0800 

0.0747 

18 

°2 

0.0802 

0.0779 

0.0739 

32 

C°2 

0.0806 

0.0804 

0.0731 

44 

NO 

0.0008 

0.0049 

0.0139 

30 

OH 

0.00006 

0.0012 

0.0078 

17 

CO 

— 

0.0003 

0.0069 

28 

O 

— 

0.0001 

0.0028 

16 

H2 

— 

— 

0.0011 

2 

'  H 

— 

— 

0.0006 

1.0 

MASS  FRACTIONS 


Species 

T  =  2.700“R 

T  =  3 . 600  °R 

T  =  4.500-R 

N2 

0.737 

0.738 

0.730 

C°2 

0.123 

0.124 

0.113 

°2 

0.089 

0.087 

0.083 

H2° 

0.050 

0.050 

0.047 

NO 

0.001 

-5 

0.005 

0.015 

OH 

0.3  x  10 

0.001 

0.005 

CO 

— 

— 

0.007 

O 

— 

— 

0.002 

H2 

— 

— 

— 

H 

— 

— 

— 
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This,  therefore,  negates  the  need  for  the  equilibrium  constants  and 
greatly  simplifies  the  numerical  computations.  However,  as  a  check 
on  this  model,  equilibrium  conditions  can  be  computed  at  different 
states,  after  the  fact,  to  confirm  this  approach.  Additional 
discussion  of  computational  procedures  for  equilibrium  chemistry  is 
presented  in  the  Appendix. 

2.1.7  Model  Problems  for  Subsonic  Combustion 

To  validate  the  chemically  reacting  program  two 
model  problems  were  accomplished.  The  first  was  the  simulation  of 
cold  flow  in  the  dump  combustor  for  a  case  in  which  experimental 
data  existed.  This  model  problem  demonstrates  the  capability  to 
simulate  the  fluid  mechanical  features.  The  second  selected  was 
the  case  for  quasi-one-dimensional  flow  with  subsonic  deflagration 
of  ten  percent  equivalence  ratio  of  premixed  methane  gas  and  air. 
This  second  case  activates  the  chemistry  features  of  the  code.  The 
results  of  these  numerical  investigations  follow. 

2. 1.7.1  Cold  Flow  Computations 

In  preparation  to  performing 

computations  for  a  chemically  reacting  flowfield,  it  is  imperative 
that  the  computer  code  have  the  capability  to  perform  computations 
for  the  non-reacting  test  case.  These  cold  flow  computations  were 
made  in  a  previous  investigation  conducted  at  UDRI  (Contract 
F33615-81-C-2078 ,  Task  019) .  Difficulty  was  encountered,  however, 
in  obtaining  agreement  with  experiment.  The  principal 
discrepancies  between  the  computed  solutions  and  the  experimental 
data  were: 

1.  The  computed  streamlines  did  not  reattach  in  the 
location  measured  in  experiment. 

2.  The  computed  average  inlet  Mach  number  was  higher 
than  the  value  calculated  in  experiment. 

3.  The  computed  frequency  of  the  self-excited 
oscillation  was  approximately  twice  the  value 
measured  in  experiment. 
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In  the  previous  investigation  attempts 
were  made  to  improve  the  accuracy  of  the  computed  solutions  by 
various  means  including: 

1.  Performing  computations  with  increased  grid 
resolutions  in  the  radial  direction. 

2.  Performing  computations  with  various  forms  of 
algebraic  turbulence  models. 

Based  on  the  results  of  the  previous 
investigation,  it  was  deemed  necessary  to  perform  a  series  of 
numerical  experiments  while  varying  the  axial  resolution  of  the 
computational  grid.  It  is  important  that  the  computational  grid 
have  sufficient  resolution  in  the  axial  directions  to  resolve  the 
large  scale  vortical  structures.  It  was  hypothesized  and  later 
verified  by  computation  that  by  increasing  the  resolution  in  the 
axial  direction  the  large  scale  vortical  structures  would  be 
resolved.  Resolution  of  the  large  scale  structures  would  allow 
more  kinetic  energy  of  the  flow  to  be  dissipated  through  vorticity. 
This  decrease  in  kinetic  energy  of  the  flow  would  result  in 
reattachment  in  a  region  closer  to  the  reattachment  point  measured 
by  experiment.  The  reattachment  point  computed  in  the  previous 
investigation  was  located  on  the  downstream  choke  nozzle  (see 
Figure  1)  (i.e.,  the  computed  field  streamlines  flow  never 
reattached  in  the  dump  combustor) . 

The  configuration  of  the  ramjet  dump 
combustor  is  shown  in  Figure  1.  This  configuration  is 
axisymmetric,  hence,  the  axisymmetric  forms  of  the  time-dependent 
Navier-Stokes  equations  can  be  used.  Note  for  the  nonreacting  case 
the  species  equations  are  dropped  from  the  analysis. 

Ut+  Ex  +  (1/r) (rF) r  =  H 
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where 


and  where  the  components  of  the  stress  tensor  are  given  by 


P 

P 

P 


The  dependent  variables  for  this  system  of  equations  are  U(p,  pu, 
pv,  pe) .  Sutherland's  viscosity  equation,  the  equation  of  state, 
and  the  Prandtl  number  are  specified  to  close  this  system  of 
equations. 
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For  the  present  investigation,  the  no¬ 
slip  boundary  conditions  are  applied  at  the  outer  duct  walls  and  a 
symmetry  boundary  condition  is  applied  at  the  duct  centerline. 

Treatment  of  the  inflow  and  outflow 
conditions  are  among  the  primary  focuses  of  this  investigation. 
These  conditions  consist  of  the  following: 

Subsonic  Inflow:  PQ,  T  ,  v  =  0  specified 

£(P- yj  _  0 
3x  “  0 

Supersonic  Outflow: 

c  ox  ox  ox  ox 

These  conditions  are  shown  in  Figure  2. 

2. 1.7. 1.1  Computational  Grid 

To  study  the  effects  of  axial  grid 
resolution  on  self-excited  oscillating  flows  a  series  of 
computations  were  performed  by  employing  grids  with  various  axial 
spacing  and  comparisons  were  made  with  experiment  for  three 
computational  grids. 

•  Grid  1-74  axial  x  41  radial  points 
(coarse) 

•  Grid  2  -  151  axial  x  41  radial  points 
(medium) 

•  Grid  3  -  301  axial  x  41  radial  points 
(fine) 

The  axial  spacing  for  the  coarse  grid  was  exponentially  stretched 
with  a  minimum  Ax  of  0.2  inch  at  the  dump  plane  and  a  maximum  Ax  of 
0.5  inch  near  the  throat  of  the  downstream  choke  nozzle.  The  axial 
spacing  for  the  medium  and  fine  grids  was  constant  at  Ax  =  0.2  and 
0.1  inch,  respectively.  All  of  the  computations  were  made  with  the 
axisymmetric  unsteady  Navier-Stokes  code  first  developed  by 

4 

Shang. 
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Figure  2.  Boundary  condition 


2. 1.7. 1.2  Numerical  Results 

Computations  were  performed  for  the 
experimental  ramjet  configuration  shown  in  Figure  1.  The  Reynolds 
number  based  on  inlet  pipe  diameter  was  1.0  x  io6.  The  average 
static  pressure  and  Mach  number  upstreams  of  the  dump  plane  were 
1,570  psf  and  0.55,  respectively.  The  total  temperature,  specified 
at  the  upstream  boundary  condition,  and  the  wall  temperature  were 
held  constant  at  522 °R. 

Computed  pressure  histograms  at 
locations  P4,  Pfi,  and  P1Q  (reference  Figure  1)  for  the  coarse, 
medium,  and  fine  computational  grids  are  shown  in  Figures  3,  4,  and 
5,  respectively.  The  results  of  the  experimental  measurements 
taken  in  the  dump  combustor  for  the  above  test  conditions  are  shown 
in  Figure  6.  The  increase  in  grid  point  density  has  a  dramatic 
impact  on  the  pressure  histograms,  especially  from  the  coarse  to 
medium  grid.  As  the  resolution  in  the  axial  direction  is  refined, 
the  computed  results  approach  the  experimental  data  (see  Table  5) 
where  A  and  f  are  the  average  peak  to  peak  amplitude  and  average 
frequency  respectively. 

TABLE  5 

SUMMARY  OF  COMPUTATIONS 


Axial  Grid  Points 

*4 

A. 

4 

f4(H2) 

M  Inlet 

71 

1,306 

125 

357 

0.77 

151 

1,405 

320 

340 

0.63 

301 

1,409 

320 

303 

0.59 

Experiment 

1,570 

202 

190 

0.55 
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Figure  3.  Pressure  hist^g 
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Figure  5.  Pressure  histograms  —  fine  grid  (301  x  U l ) 


Another  benefit  of  the  increased 
resolution  of  the  medium  and  fine  grids  is  the  ability  to  better 
observe  the  shedding  of  vorticies  generated  in  the  dump  region. 

The  shedding  phenomena  is  demonstrated  in  Figure  7.  The  time 
history  of  the  vorticity  clearly  depicts  the  unsteady  character  of 
the  flow  and  also  provides  new  insight  into  the  flow  reattachment 
phenomena . 

In  the  previous  investigation  it  was 
noted  that  for  the  laminar  test  cases  the  computed  flowfield  did 
not  reattach  in  the  combustor  region.  The  experimental  data 
indicated  an  average  reattachment  point  near  x  *  11.0  inches. 

When  the  computations  were  done  on  the  medium  and  fine  grids,  the 
reattachment  location  was  in  agreement  with  experiment.  Figure  8 
displays  how  the  reattachment  points  move  with  time.  Note  that 
there  are  several  reattachment  points  at  any  one  time,  but  the 
primary  reattachment  point  varies  from  8.0  inches  to  approximately 
13.5  inches. 

The  results  of  the  numerical  study 
were  very  encouraging.  The  axial  resolution  of  the  computational 
grid  has  a  dramatic  impact  on  the  accuracy  of  unsteady  flow 
computations.  It  has  been  demonstrated  that  increasing  the  grid 
point  density  in  the  axial  direction  leads  to  better  resolution  of 
the  large  scale  vortical  structures.  This,  in  turn,  leads  to 
increased  dissipation  and  the  computed  results  more  closely  agree 
with  experiment.  However,  the  accuracy  of  the  computations  still 
could  be  further  improved. 

Over  the  course  of  the  numerical 

investigation  it  became  apparent  that  in  order  to  obtain  reasonable 
results  the  computer  code  has  to  be  capable  of  adequately  modeling 
the  viscous  dissipation.  The  accuracy  of  the  computed  solutions 
was  greatly  increased  when  the  grid  was  refined  to  allow  for 
resolution  of  the  large  scale  eddies.  It  is  conjectured  that  the 
accuracy  of  the  computer  code  could  be  further  increased  by 
modeling  the  fine  scale  turbulence  in  the  shear  layer  by 
application  of  a  "law  of  the  wake”  turbulence  model.  Another 
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Figure  8.  Reat tachment  phenomenon. 


important  consideration  is  the  three-dimensional  effect.  By 
performing  the  computations  in  the  three-dimensions  instead  of 
assuming  an  axisymmetric  flow,  the  three-dimensional  effect  of 
turbulence  may  be  better  simulated. 

These  results  suggest  guidelines  for 
the  minimum  number  of  grid  points  necessary  for  resolving  the  time- 
dependent  flow  features. 

2. 1.7. 2  Ouasi-One-Dimensional  Model  Problem 

The  quasi-one-dimensional  model  problem 
consisted  of  subsonic  deflagration  of  methane  gas  and  air.  The 
methane-air  mixture  was  heated  above  the  ignition  point  of  the 
mixture  by  passing  through  a  normal  shock.  The  numerical  model 
included  additional  air  so  the  computational  fluid  represented  a 
10  percent  equivalence  ratio  of  the  gas  mixture. 

2. 1.7. 2.1  Governing  Equations 

The  chemical  formula  for  the 
stoichiometric  combustion  of  methane  in  air  is: 

Kf 

CH.  +  20o  +  7.52  N0  -  CO_  +  2H_0  +  7.52  N_ 

4  2  2  2  2  2 


The  governing  fluid  dynamic  equations  are: 


PA 

puA 

puA 

+ 

{pu2  +  P  -  Txx}A 

pEtA 

{u(pEt  +  P  -  Txx)  +  qx}A 

dY. 

PY.A 

t 

<'<UYi  -  “dT  »A 

{  p  -  j  }  — 

1  *  xxr  dx 
1  xx  ^xJ  dx 

’  a  .ndYi  dA 

mi A  -  pD^r  d^ 
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where 


F  =  h  +  _  £ 

Et  h  2  p 


and 


T 

h  =  I  vj  J 

Tref 


CP.jdT 


h°) 


The  primary  difficulties  encountered 

in  the  integration  of  Equation  (2)  is  the  representation  of  the 

mass  production  terms  (m^) .  If  a  global  one-step  model  is  used  to 

calculate  the  rate  of  fuel  consumption,  all  other  mass  production 

rates  are  algebraically  determined  via  the  law  of  mass  action. 

5 

Coffee,  et  al.,  developed  an  expression  for  fuel  consumption. 

They  proposed: 


m 


fuel 


Yfuel 


^  Yoxidizer^ 


where  i/p  and  vq  are  the  stoichiometric  coefficients  for  the  global 
reaction  formula.  The  constants  A.  and  T  are  determined  from  best 

-L  d 

fit  curves  generated  by  a  detailed  reaction  rate  model  which 
employs  experimental  data  as  input. 

2. 1.7. 2. 2  Initial  and  Boundary  Conditions 

The  incoming  flow  was  supersonic, 

therefore,  all  flow  variables  were  specified  at  station  (1) ,  Figure 
9.  Since  the  outflow  was  subsonic,  the  back  pressure  | P2 |  was 
specified  to  ensure  the  formation  of  a  normal  shock  in  the  nozzle. 

The  incoming  flow  contained  10  percent 
of  the  stoichiometric  mass  fraction  of  methane.  By  a  simple 
analysis,  it  was  determined  that  only  13  percent  of  the 
stoichiometric  mass  fraction  could  be  burned  before  thermal  choking 
would  occur. 
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mx  =  l.a 

Pj  =  1000.0  Psf. 
Tj  =  1350.0  R 
Y1  =  Constant 


Figure  9.  Shock  Induced  Deflagration  Model  Problem. 
Quasi-One-Dimensional . 
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Since  nitrogen  was  the  dominant 
species,  the  fluid  properties  were  assumed  to  be  those  of 
nitrogen. 

=  Constant  =  5.830  x  io3  ^ f£, 

P  slug  -  R 


k 


2.875  X  10 


-3 


0.94 


4  9 1  °  RJ 


Lbf  -  ft 
ft  -  S  -  R* 


D 


n2,ch4 


=  2.1527  x  10 


-4 


fti 


A 


1 


-1.0  x 


1018  {- 


ft' 


slug" 


1 

S 


2. 1.7. 2. 3  Results 

Figures  10  and  11  depict  axial 

pressure  through  the  nozzle  for  the  reacting  and  nonreacting  cases. 
The  test  case  shown  in  Figure  10  incorporated  larger  amounts  of 
normal  stress  damping6  than  the  nonreacting  calculation.  The  test 
case  shown  in  Figure  12  has  the  same  amount  of  normal  stress 
damping  as  the  nonreacting  calculation.  From  the  above,  it  was 
observed  that  normal  stress  damping  adequately  removes  shock 
induced  oscillations  without  affecting  the  jump  conditions. 

Figures  13  and  14  depict  axial  temperature  distribution  for  the 
reacting  and  nonreacting  cases.  The  species  mass  fractions  for  the 
reacting  case  are  shown  in  Figure  15. 

The  addition  of  the  chemical  rate 

equations  did  not  reduce  the  stability  of  the  MacCormack  algorithm. 
The  time  step  used  was  80  percent  of  the  inviscid  CFL  limit.  This 
behavior  is  opposite  to  that  reported  by  Eklund,  et  al.,7  and  was 
primarily  due  to  the  global  reaction  rate  model. 

2 . 2  Supersonic  Combustion 

The  second  major  effort  in  the  computation  of  ramjet 
internal  flowfields  was  the  investigation  of  supersonic  combustion. 
For  this  supersonic  flow  regime  advantage  may  be  taken  of  an 
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AXIAL  PRESSURE  —  10%  STOIC. 


AXIAL  POSITION  —  FEET 


AXIAL  PRESSURE  —  COLD  FLOW 
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AXIAL  POSITION  ---  FEET 

Figure  11.  Axial  Pressure  Profile  (Nonreaeting  Flow) 
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MASS  FRACTIONS  —  10%  STOIC. 
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AXIAL  POSITION  —  FEET 

il  Mass  Fractions  (Reacting  Flow). 


approximation  which  cannot  be  utilized  in  subsonic  unsteady  flows, 
i.e.,  PNS.  The  parabolized  Navier-Stokes  equations  are  appropriate 
for  only  steady,  supersonic,  attached  flows.  This  limitation, 
however,  results  in  at  least  an  order  of  magnitude  reduction  in 
computer  time.  The  savings  in  computer  resources  is  extremely 
desirable  in  systems  studies  for  the  Aerospace  plane.  Therefore, 
supersonic  combustion  of  hydrogen  fuel  was  examined  using  the  PNS 
method . 

In  the  following  sections  the  governing  equations, 
boundary  conditions,  and  chemistry  for  hydrogen/air  reaction  will 
be  presented.  Model  problems  are  solved  for  several  cases  to 
validate  the  method. 


2.2.1 


Governing  Equations 


The  steady  parabolized  Navier-Stokes  equations 
were  used  in  this  investigation.  In  divergence  form,  the  steady 
Navier-Stokes  equations  are: 


E  +  F  +  G  =  H 
x  y  z 


Where  the  E,  F,  G,  and  H  vectors  are  defined  as: 


E  = 


fiU 

2 

p  u  +  p  -  r 


xx 


p  UV  -  T 


xy 


puw  -  Txz 
(Et  +  p)U  -  UT 
puYj  -  „DYiiJ£ 


-  vr  -  wr  -  q 
xx  xy  xz 
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pv  +  p  -  T, 

PVW  -  T 


(Et  +  P)V  -  UTxy 
pVYi  -  />DY.>y 


puw  -  T 


PVW  -  T 


/,W  +  p  -  rz2 

(Et  +  p)W  -  UT, 

„wy.  -  ,mitt 


The  governing  equations  were  then  transformed 
from  the  physical  to  the  computational  plane  to  simplify  the 
numerical  integration  algorithm.  The  transformation  used  in  this 
study  was  as  follows: 

f  =  f(x) 
n  -  V (x,y, z) 
f  ■  f (x,y,z) 


Applying  the  above  transformation  to  the  vector 


equation  yields: 


+  r?  E  +  n  F  +  f ,.F_  +  r?  G  +  CGr)  i  =  9 


where  J  is  the  Jacobian  of  the  transformation.  The  equations  are 

then  parabolized  in  the  transformed  plane  to  insure  maximum 

accuracy.  In  addition,  the  splitting  of  the  axial  streamwise 

pressure  gradient  was  performed  as  specified  by  Vigneron,  Rakich, 
8 

and  Tannehill.  Performing  the  parabolization  and  recasting  the 
above  equation  into  divergence  form  results  in: 


Where  the  E*  and  the  P  vector  and  the  scaler  u  are  defined  as 
follows  (Vigneron,  et  al.): 


E*  =  pu 

2 

pu  +  up 

puv 

puw 

u(Et  +  p) 
upYi 

P  =  0 

(l-w)p 

0 

0 

0 


% 


w 


(l+(7-l)M2) 


The  above  equations  contain  the  geometric 

9 

conservation  law  (GCL)  terms  (Thomas  and  Lombard  )  and  the  metric 
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derivatives  which  result  from  returning  to  the  conservation  law 
form  in  the  transformed  plane.  It  was  found  that  the  GCL  terms 
must  be  included  into  the  formulation  in  order  to  eliminate  source 
term-like  errors  (Gielda  and  McRae10'11). 

Applying  the  transformation  and  parabolization  to 
the  viscous  stress  and  heat  flux  terms  yields  the  following 
relations: 

rxx  “  3*  <2<Vx  +  V*’  ■  <Vy  +  Vy  +  Vz  +  Vz)l 

ryy  ■  “  {2<Vy  +  Vy>  '  <Vx  +  Vx  +  V*  +  Vz)} 

TZZ  =  3^  <2<Vz  +  Wffz)  "  (Vn^y  +  Vrfy  +  U„1X  +  UffX>t 

TXy  ”  l‘^un’ly  +  uffy  +  Vx  +  Vx> 

Tyz  -  ^Vy  +  Vy  +  V*  +  V*> 

rxz  -  "<Vz  +  V*  +  Vx  +  Vx> 


^x  -  <Vx  +  Tffx>  +  'D  Jhi<Yi,„"x  +  Yi,ffx' 


^y  =  Pr  y  +  Vx>  +  pD  2hi(Yi,r?,7x  +  Yi,ffx) 


5z  “  -fe  +  Vz}  +  0°  Zhi<Yi,„”z  +  Yi,ffz> 


These  relations  were  then  used  in  the  formulation 
of  the  E,  F,  and  G  vectors.  The  thermodynamic  properties  (Cp,  R, 
etc.)  were  functions  of  species  mass  fraction  and  the  local  static 
temperature.  Cp  was  determined  from  the  local  mass  fraction  and 
Cpi  were  calculated  via  curve  fits  with  ranges  of  applicability  up 
to  6500*R. 
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2.2.2 


Numerical  Method 


12 

MacCormack's  well-known  predictor-corrector 
scheme  was  used  to  integrate  the  governing  equations.  The 
application  of  this  numerical  integration  scheme  is  well-documented 
in  References  10  and  11,  and,  therefore,  will  not  be  repeated. 

The  metrics  in  the  previous  equations  were 
calculated  by  second  order  central  differences  except  for  the  f 
derivatives  which  were  calculated  by  first  order  backward 
differences.  First  order  backward  differences  for  the  f 
derivatives  were  required  since  only  two  planes  of  data  are  stored 
in  computer  memory  for  each  marching  step.  After  each  marching 
step  the  current  information  at  the  n+1  level  is  shifted  to  the  nth 
level . 

The  viscous  and  heat  flux  terms  used  in  forming 
the  E,  F,  and  G  vectors  were  differenced  in  the  opposite  direction 
from  the  flux  derivative.  For  example,  in  forming 

{E77x  +  F  rjy  +  cr)2}  1/J 

in  the  predictor  step,  the  derivatives  with  respect  to  rj  were 
calculated  by  one-sided  backward  differences,  while  central 
differencing  was  used  for  the  f  derivatives.  This  must  be  done  to 
maintain  the  three-point  stencil  for  the  viscous  terms.  Likewise, 
in  forming 

{Efx  +  Ffy  +  Gfz}  1/J 

in  the  predictor  step,  the  f  derivatives  were  calculated  via  one¬ 
sided  backward  differences  and  the  r]  derivatives  were  calculated  by 
central  differences.  In  the  corrector  sweep  the  one-sided 
differences  switch  from  backward  to  forward. 

Streamwise  Pressure  Gradient 

As  stated  in  the  section  on  governing  equations,  the  streamwise 
pressure  gradient  is  split  in  a  manner  identical  to  that  proposed 
by  Vigneron,  et  al.  Vigneron  has  proved  that  a  fraction  of  the 
streamwise  pressure  gradient  can  be  included  in  the  governing 
equations  and  still  maintain  algorithm  stability  while  marching 
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through  regions  of  locally  subsonic  flow.  Recalling  the  pressure 
splitting  scheme,  where  E*  and  the  P  vector  are  defined  as: 


pu 

p  = 

0 

pu2  +  up 

(l-w)p 

puv 

0 

puw 

0 

(Et  +  P)U 

0 

PUY. 

0 

The  function  u  is  the  maximum  allowable  fraction  of  the  streamwise 
pressure  gradient  which  can  be  retained  and  maintain  stable 
marching.  Taking  into  account  the  nonlinear  effects  of  the 
governing  equations,  u  is  redefined  as: 

a-yM2 

u  =  - - — 

i+cr-i)M2x 

The  value  of  the  safety  factor  o  was  approximately  0.75. 

The  major  benefit  of  the  Vigneron  technique  is  its  effect  on 
the  decoding  of  the  axial  flux  vector  E*.  A  quadratic  equation  for 
u  velocity  component  results  when  the  dependent  variables  are 
computed  from  the  components  of  the  E*  or  E*  vectors.  The  roots  of 
the  quadratic  equation  correspond  to  either  the  subsonic  or 
supersonic  solution  (elliptic  or  hyperbolic  character  of  the 
governing  equation) .  By  incorporating  the  Vigneron  scheme,  the 
character  of  the  governing  equations  remains  hyperbolic  throughout 
the  flowfield.  Therefore,  only  the  root  corresponding  to  the 
hyperbolic  solution  is  used. 

2.2.3  Initial  and  Boundary  Conditions 

The  initialization  of  a  PNS  code  requires  an 
initial  starting  solution.  The  PNS  solution  is  then  marched 
downstream  from  this  initial  data  plane.  During  the  investigation 
two  types  of  initial  conditions  were  employed:  (1)  The  initial 
data  plane  was  set  equal  to  free  stream  conditions  or  the 
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conditions  upstream  of  the  inlet.  (2)  The  initial  data  plane  was 
generated  by  another  integration  scheme  (i.e.,  a  time  dependent  FNS 
solution) .  The  latter  technique  is  most  beneficial  when  regions  of 
streamwise  separation  are  encountered. 

At  all  solid  surfaces,  isothermal  and  no-slip 
boundary  conditions  were  applied.  The  pressure  and  species  mass 
fractions  were  determined  from  zero  dp/dn  and  zero  dY./dn  boundary 
conditions  respectively. 

2.2.4  Numerical  Smoothing 

The  MacCormack  explicit  algorithm  is  second  order 

accurate  in  the  marching  as  well  as  the  transverse  and  normal 

coordinate  direction.  Since  the  algorithm  is  second  order 

accurate,  the  leading  term  is  third  order  or  dispersive.  The 

MacCormack  explicit  PNS  code  incorporated  both  normal  stress 
6  13 

damping  (McRae  )  and  modified  MacCormack  second  order  smoothing 
(Gielda  and  McRae10'11).  Figures  16  and  17  depict  the  surface 
pressure  on  a  secant-ogive-cylinder  projectile  at  Mach  4  and  angle 
of  attack,  without  and  with  normal  stress  damping.  The 
oscillations  associated  with  the  program  start-up  are  observed  to 
be  suppressed  immediately  with  no  offset  in  surface  pressure.  This 
form  of  damping  results  when  the  normal  stress  terms  are  modified 
as  shown  below. 


defined  as: 


The  normal  stress  terms  (txx,  ryy,  and  rzz)  are 


r 


xx 


=  0(2ux) 


+  % 


where  a is  normally  defined  as: 

~  §M(V*V) 

Normal  stress  damping  results  when  a is  multiplied  by  a  factor  P, 
where  P  is  of  the  form: 

p 

P  =  ( 1-DNS1  -  DNS 2  5~) 

Poo  Poo 
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u  *p'rri  } . \  i"1 1"  i  »"■  \  i  '[  '[  t  "v  rj-n  \  i  {  i  i  \  >  [ 

0.04  0.05  0.06  0.07  0.08  0.09  O.LQ 

X/L  REL  =  6.02E+06 


Figure  16.  Computed  surface  pressure  on  a  secant-ogive-cylinder- 
boat-tail  projectile  at  4.0  decrees  angle  of  attack 
without  normal  stress  damping  {(^*4.0). 
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PU/PINF 


Triangle  —  Windward  Plana. 
Diamond  —  Mid  Plana. 
Square  —  leeward  Plana. 


X/L  DMS1  =  1 


DMS2  =  100 


Figure  17.  Computed  surface  pressure  on  a  secant-ogive-cylinder- 
boat-tail  projectile  at  4.0  degrees  angle  of  attack 
with  normal  stress  damping  (M  *4.0). 


46 


The  typical  values  of  DNSl  and  DNS2  varied  from  0<DNS1<10  and 
0<DNS2<10 . 

In  addition  to  normal  stress  damping,  second 
order  smoothing  was  also  used.  The  functional  form  of  this  damping 
was: 

Damping  =  (A1  ( j  ,k)  E(  j  ,k,n)  +  A2 ( j , k) E ( j , k, n) )  (-J) 

Where  A1  and  A2  are  defined  by: 

Al(j,K)  -  CDY  ABS  {p  ,  j  +  1 ,  ^ !p +p  ,  j -1,  k,  > 

A2  ( j  ,k)  =  CDS  ABS  {p(jiit+ij-;-2p^k)4-p(j  ,k-l)} 

The  second  order  explicit  damping  was  applied  in 
both  the  predictor  and  corrector  steps.  In  the  predictor  the 
damping  was  based  on  E*n,  while  in  the  corrector  the  damping  was 
based  on  E*n+1.  The  values  of  CDY  and  CDZ  varied  from  0.0  to  0.05. 
The  second  order  damping  was  used  sparingly  and  primarily  during 
program  initiation. 

Once  the  initial  starting  transients  were 
eliminated  the  normal  stress  and  second  order  damping  were  removed. 
Typically  this  was  accomplished  by  reducing  the  coefficients  DNSl, 
DNS 2 ,  CDY,  and  CDZ  by  0.5  percent  per  marching  step. 

2.2.5  Reaction  Rate  Modeling 

Hydrogen-oxygen  reaction  has  been  thoroughly 

investigated.  Twelve  species  and  25  reactions  have  been  identified 

14 

and  are  listed  in  Table  6.  A  simpler  8-reaction  system  is  also 
identified  in  the  Table.  However,  to  begin  the  investigation  a 
single  global  reaction  was  explored. 

For  the  preliminary  supersonic/hypersonic 
combustion  computations  global  chemical  reaction  rates  were 
employed.  The  global  reaction  rate  models  were  incorporated  in  an 
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effort  to  reduce  computational  and  storage  requirements.  Even 
though  the  global  models  do  not  predict  accurate  flame  structure, 
they  do  perform  adequately  in  computing  the  jump  conditions  across 
the  flame.  The  global  reaction  for  the  stoichiometric  combustion 
of  hydrogen  in  air  is  given  by: 

2H2  +  02  +  3.76N2 - -  2H20  +  3.76N2 

The  global  reaction  rate  Varma,  Chatwani,  and  Bracco  was  used  in 
this  investigation.  In  addition  to  the  finite  rate  model,  an 
additional  reaction  rate  model  based  on  chemical  equilibrium 
(Varma,  et  al.)  was  considered. 

The  equilibrium  reaction  rate  models  are  well- 
suited  for  high  temperature  flames  encountered  in  hydrogen 
combustion.  These  models  allow  for  forward  and  backward  reaction 
rates.  Unfortunately,  the  equilibrium  models  are  not  well-suited 
for  vector  processing  computers  since  they  require  frequent  table 
"look-up"  operations. 

2.2.6  Finite  Rate  Model 


As  stated  above,  the  global  reaction  rate  model 
of  Varma,  et  al.,  was  considered.  The  global  finite  rate  equation 
is  written  as: 


m 


fuel 


=  -ax[h2 


VH2> 


The  quantities  in  the  brackets  represent  the  local  concentration  of 

the  species.  Mw(H2)  is  the  molecular  weight  of  hydrogen.  The 

activation  temperature,  T  ,  and  the  pre-exponential  constant 

were  found  to  be  strong  functions  of  the  equivalence  ratio  of  the 
15 

mixture. 

2.2.7  Equilibrium  Reaction  Rate  Modeling 

The  primary  limitation  of  the  finite  rate 
chemistry  models  are  that  they  are  defined  by  a  single  Arrhenius 
type  reaction.  However,  in  regions  of  high  temperature  these 
models  do  not  perform  as  well  since  there  are  significant  shifts  in 
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TABLE  6  SPECIES  ANO  REACT  IONS* 


Species 

1 

H" 

7 

N2 

2 

0 

8 

N 

3 

h2o 

9 

NO 

4 

OH 

10 

no2 

S 

°2 

11 

ho2 

6 

H2 

12 

hno2 

Reactions 


1 

HHOj  +  H  -  NO  ♦  OH  ♦  H 

14 

OH  ♦  OH  -  H  +  H02 

2 

N02  ♦  H  -  NO  ♦  0  +  H 

15 

H20  +  0  -  H  +  H02 

*3 

h2*h-h  +  h*h 

16 

OH  ♦  02  -  0  ♦  H0Z 

*4 

o2*h-o  +  o  +  h 

17 

H20  +  02  -  OH  ♦  H02 

*5 

H20  ♦  H  -  OH  ♦  H  ♦  H 

18 

H20  +  OH  -  H2  ♦  H02 

*6 

OH+H-O+H+H 

19 

0  +  N2  -  N  +  NO 

7 

H02  ♦  H  -  H  ♦  02  +  H 

20 

H  ♦  NO  -  N  +  OH 

*8 

H20  +  0  -  OH  ♦  OH 

21 

0  +  NO  -  N  ♦  02 

*9 

H20  ♦  H  -  OH  ♦  H2 

22 

NO  +  OH  -  H  +  N02 

*10 

02  ♦  H  -  OH  ♦  0 

23 

NO  +  02  -  0  +  N02 

*11 

H?  ♦  0  -  OH  +  H 

24 

N02  ♦  H2  -  H  ♦  HN02 

12 

H2  ♦  Oj  -  OH  ♦  OH 

25 

N02  +  OH  -  NO  +  H02 

13 

H2  +  02  -  H  +  H02 

Reaction  Rates 


Forward  Rate  Constant**  Rever  ;e  Rate  Constant1* 


A 

8 

c 

A 

8 

c 

1 

5.0 

X 

rstoflOCO- icoiown^nnn^-^- •  — 

—W  —+  OJ  1  •— «  —»  ~  «— •  *-*•  •  •— «  •— «  •— * 

ooooooooooooooooooooooooo 

-1.0 

25000 

8.0 

X 

0 

-1000 

2 

1.1 

X 

0 

32712 

1.1 

X 

0 

-941 

3 

5.5 

X 

-1.0 

51987 

1.8 

X 

htfm 

0 

4 

7.2 

X 

-1.0 

59340 

4.0 

X 

0 

5 

5.2 

X 

-1.5 

S9386 

4.4 

X 

iiStltfli 

0 

6 

8.5 

X 

-1.0 

50830 

7.1 

X 

IBPIV 

0 

7 

1.7 

X 

0 

23100 

1.1 

X 

mKVSm 

-440 

8 

5.8 

X 

0 

9059 

5.3 

X 

503 

9 

8.4 

X 

0 

10116 

2.0 

X 

2600 

10 

2  2 

X 

8455 

1.5 

X 

0 

11 

7.5 

X 

■■■ 

5586 

3.0 

X 

4429 

12 

1.7 

X 

24232 

5.7 

X 

14922 

13 

1.9 

X 

BIS 

24100 

1.3 

X 

0 

C 

14 

1.7 

X 

BOTH 

21137 

6.0 

X 

0 

0 

IS 

5.8 

X 

iHWi 

28686 

3.0 

X 

0 

0 

16 

3.7 

X 

27840 

1.0 

X 

0 

0 

17 

2.0 

X 

■m 

36296 

1.2 

X 

0 

0 

18 

1.2 

X 

0.21 

39815 

1.7 

X 

0 

12582 

19 

5.0 

X 

0 

37940 

1.1 

X 

0 

0 

20 

1.7 

X 

mam 

24500 

4.5 

X 

0 

0 

21 

2.4 

X 

■n 

19200 

1.0 

X 

0.5 

3120 

22 

2.0 

X 

15500 

3.5 

X 

740 

23 

1.0 

X 

0 

22800 

1.0 

X 

iBfl. 

302 

24 

2.4 

X 

0 

14500 

5.0 

X 

m£m 

1500 

25 

1.0 

X 

0.5 

6000 

3.0 

X 

0.5 

1200 

*The  first  7  species  and  the  starred  reactions  constitute  the  8-reaction 
system. 

k.  a  s 

For*  of  rate  constant  It  k  ■  AT°exp(-C/T)  with  It  in  cm  /mole-sec  or 
cia®/*ole2-sec. 
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chemical  equilibrium.  In  order  to  improve  on  the  accuracy  of  the 
finite  rate  models,  the  equilibrium  reaction  rate  of  Varma  was 
incorporated  into  a  version  of  the  PNS  code  (SSCPNSEQ)  modified  to 
perform  calculations  with  equilibrium  chemistry  in  a  zone  of 
supersonic  combustion. 

The  global  chemical  reaction  equation  employing 
the  equilibrium  reaction  rate  model  is  written  as: 

_ -j _ 

2H2  +  02  ♦  3.76N2  ^  2H20  +  3.76«2 


The  rate  of  production/depletion  of  species  j  is  defined  by: 

[Cj]  =  {[Cj]  -  [ C*]}/T 

where  [C^]  is  the  local  concentration  of  species  j  and  [C*^]  is  the 
concentration  of  species  j  if  the  mixture  were  in  chemical 
equilibrium.  r  is  defined  as  the  characteristic  time  of  the  global 
reaction  and  is  given  by: 

1/7  =  B[H2]0,6[O2]0*6  e(_Ta/T) 

As  in  the  case  of  the  finite  rate  models  the  constants  T  and  B  are 

a 

functions  of  the  equivalence  ratio  of  the  mixture.  For  complete 
details  see  Reference  15. 

2.2.8  Equilibrium  Chemistry 

To  devise  a  combustion  model  it  is  necessary  to 
examine  the  equilibrium  state  at  the  final  condition.  To 
accomplish  this  task  nine  species  are  included.  At  the  temperature 
levels  involved  3  species  may  be  eliminated  from  the  12  variables 
listed  in  Table  6  (i.e.,  NC>2,  H02,  and  HN02). 

The  procedure  is  similar  to  the  approach  taken  in 
Section  2.2.6.  The  mole  fractions  of  the  equilibrium  composition 
of  the  products  of  combustion  of  H2  are  determined  using  the 
equations  involving  the  equilibrium  constants  with  the  assumption 
that  the  component  gases  obey  the  equation  of  state  for  an  ideal 
gas . 
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Pi  RT 

p.  =  — * — 
i  M . 


The  general  chemical  reaction  of  the  nine  species  may  be  written  as 
follows: 

HJ  +  V  <°2  +  3'76N2>  -  nN2<N2>  +  n02(°2>  +  nH20<H20) 

+  n^NO)  +  nOH(OH)  +  nH(H)  +  n0(0)  +  nN(N)  +  nR  (H2) 

where  n^  «  the  number  of  moles  of  product  "i" 

0  =  equivalence  ratio 

Three  elements  are  involved  (H,0,N)  and,  hence, 
three  conservation  of  chemical  species  equations  may  be  specified. 


H!  2  =  \o  +  "°H  +  "H  2"H2 


a  ]  2nn  +  n„  +  nM  +  n_„  +  n . 


02  “H20  NO  ‘OH  0 


N:  3.76  *=  2nN^  +  nN0  +  nN 


Since  there  are  nine  unknown  molar  species  values,  six  equilibrium 
constants  are  required  to  close  the  system. 
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where  K  .  =  f(t)  from  JANNAF  Tables.  Two  additional  variables 
P  f  1 

appear  in  these  equations,  i.e.,  p  and  n. 
p  =  2p^  =  specified  pressure 

n  =  2n^  =  total  number  of  moles  of  the  combustion  products 
which  becomes  an  auxiliary  equation. 

The  problem  is  now  formulated.  Given  $ ,  p,  and  T 
there  exist  10  algebraic  equations  involving  9  molar  values  (n^) 
plus  the  value  for  the  total  number  of  moles  (n) .  Conceptually, 
this  system  of  equations  could  be  reduced  to  a  single  high  degree 
polynomial.  Physical  reasoning  requires  that  only  one  positive 
real  root  exist.  However,  in  place  of  reducing  the  system  to  one 
single  equation,  in  practice  the  entire  system  of  equations  is 
solved  by  iteration. 

An  example  problem  has  been  solved  for  <f>=l  and 
p=l  atmosphere,  over  the  temperature  range  of  interest.  These 
conditions  are  typical  of  the  scram jet  combustor  state  and  serve  to 
identify  the  dominant  species.  A  summary  of  the  results  are  shown 
in  the  following  table. 
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Once  the  equilibrium  state  has  been  identified 
the  overall  heat  release  during  combustion  may  be  computed  from 
the  heats  of  formation  of  the  individual  species. 

Ah  =  Sn.AHl 
c  it 

The  following  heat  of  combustion  values  were  determined  using 
information  from  Tables  7  and  8  and  JANNAF  values  for  AH£. 


T  (‘R) 

Ah  (kcal/mole) 

3600 

-59.5 

4500 

-53.7 

5400 

-27.4 

6300 

+  40.6 

From  these  results  it  is  apparent  that 
combustion  temperatures  must  be  limited  to  values  below  5000 °R  in 
order  to  maintain  useful  values  for  the  heat  of  combustion 
(Figure  18) . 

2.2.10  Model  Problems  for  Supersonic  Combustion 

During  this  investigation  two  model  problems 
were  computed.  In  order  to  accomplish  this  task  the  computed 
solutions  were  compared  to  experimental  data  and  the  computations 
of  others. 

The  validation  test  cases  consisted  of:  (1) 
Laminar,  two-dimensional  hypersonic  flow  over  a  15  degree 
compression  ramp.  (2)  Turbulent,  two-dimensional  supersonic 
combustion  of  hydrogen  in  vitiated  air. 
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TABLE  7 

HYDROGEN-AIR  COMBUSTION 
0=1  Equilibrium  States  p=l  Atmosphere 


T 

3600  ° R 

4500  °R 

5400  0  R 

6300 “R 

X0 

3467  Vf 

167  vT 

22.0  Vf 

5.15  Vf 

0.020 

0.0593 

0.122 

0.204 

“n 

0 

0 

0 

2.2lxlO-4V® 

V  n 

ko 

6.64xl0-4^ 

fn 

p 

0.0144 

v  P 

0-112 Vf 

0.490 

v  P 

koh 

0.581 

0.899 

1.20 

1.42 

“h 

1.62xl0~3-y 

nr 

p 

0.0251V? 

P 

0.1577/°" 

0.5877/f 

nN2 

1.879 

1.875 

1.860 

1.840 

nH 

0 

0.011 

0.129 

0.639 

"h2 

0.009 

0.066 

0.214 

0.319 

no 

0 

0.003 

0.049 

0.274 

nH20 

0.990 

0.912 

0.653 

0.248 

nOH 

0.003 

0.033 

0.136 

0.233 

"°2 

0.003 

0.020 

0.  060 

0.084 

nNO 

0.001 

0.011 

0.041 

0.080 

-4 

nN 

0 

0 

0 

6.10 

n 

2.866 

2.931 

3 . 142 

3.717 
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TABLE  8 

HEATS  OF  FORMATION 


Species 

ni 

AH£  kcal/mole 

ni 

2ru  AH£  kcal/mole 

3600 

H 

0 

— 

0 

1 

0 

0 

— 

0 

H2° 

0.990 

-60.15 

-59.5 

OH 

0.003 

8.8 

0.03 

i 

NO 

0.001 

21.6 

0.02 

-59.5 

i 

4500 

H 

0.011 

54.6 

0.6 

O 

0.003 

61.2 

0.2 

h20 

0.912 

-60.36 

-55.0 

OH 

0.033 

8.6 

0.3 

-53.7 

NO 

0.011 

21.6 

0.2 

5400 

H 

0.129 

54.9 

7.08 

O 

0.049 

61.4 

3.01 

H2° 

0.653 

-60.53 

-39.53 

-27.4 

OH 

0.136 

8.4 

1.14 

NO 

0.041 

21.5 

0.88 

! 

6300 

H 

0.639 

55.16 

35.25 

0 

0.274 

61.45 

16.84 

H2° 

0.248 

-60.70 

-15.05 

OH 

0.233 

8.17 

1.90 

+  40.6 

NO 

0.080 

21.33 

1.71 
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HYDROGEN -AIR 


Cl  "‘bust  ion 


2.2.10.1 


a  =  0.75  CDY  =  0.025 

The  computational  grid  consisted  of 
51  normal  grid  points  geometrically  spaced  from  0.0  to  0.129L. 

As  the  computation  was  marched  up  the  compression  ramp,  the  grid 
was  shifted  by  Arj  where 

Arj  =  Ax  Tan(0w) 

The  first  point  away  from  the  body 
surface  was  located  at  a  distance  of  Ay1/L=0.002  and  the  marching 
step  size  ranged  from  one-quarter  to  three-quarters  Ay^L.  The 
SSCPNS  code  has  proven  to  be  extremely  efficient  in  computing 
hypersonic  flowfields.  The  above  computation  was  performed  in 
2.2  seconds  on  a  CRAY  X-MP  computer.  This  compares  to  a 
computation  time  of  1.72  seconds  for  the  implicit  calculation  of 
Lawrence . 
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n  =  14.1 

00 


Figure  19.  Ramp  conflgurati 


For  this  test  case,  normal  stress 
damping  was  introduced  at  solution  initialization  and  then 
reduced  by  0.5  percent  per  step.  In  addition,  a  small  amount  of 
second  order  smoothing  was  introduced  to  reduce  oscillations. 

The  following  damping  coefficients  were  used: 

DNS1  =  DNS 2  =0.25 


Figure  20  depicts  the  surface  pressure  coefficient  on  the  ramp 
surface,  as  a  function  of  x,  where  is  defined  as: 


cp- 


P  V 
r  oo  oo 


Good  agreement  was  obtained  with  Holden's  experimental  data  and 
the  FNS  and  PNS  computations  of  Hung  and  MacCormack  and  Lawrence, 
et  al.,  respectively.  The  Stanton  number  is  plotted  as  a 
function  of  x  in  Figure  21;  again  note  good  agreement  was 
obtained. 

2.2.10.2  Experiment  of  Burrows  and  Kurkov 

In  the  second  test  case  the  SSCPNS 
code  was  used  to  compute  turbulent,  two-dimensional,  supersonic 

combustion  flowfields.  The  test  case  conditions  are  those  of  the 

.  .  .  .  .  19 

supersonic  combustion  investigation  of  Burrows  and  Kurkov.  A 

schematic  of  their  experiment  is  shown  in  Figure  22.  For  their 

investigation,  Mach  2.44  vitiated  air  is  combusted  with  hydrogen. 

The  hydrogen  was  injected  with  sonic  speed  through  the  small 

injection  slot. 

The  computational  grid  used  for  this 

test  case  consisted  of  81  points  in  the  normal  direction.  The 

grid  spacing  in  the  normal  direction  was  exponentially  stretched 

to  insure  adequate  resolution  of  the  shear  layer.  A  simple 

2  0 

algebraic  turbulence  model  (Baldwm-Lomax  )  was  incorporated  tc 
model  the  turbulent  mixing. 
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LAWRENCE'S  CALCULATION 


function  of  X.  Holden's  test  conditions 


STANTON  NO. 


•  LAWRENCE'S  CALCULATION 

- HAC  PNS 

°  hunTmo  mIccomuck  NEAT  TRANSFER  COEFF. 


X/L 


Figure  21.  Comparison  of  heat  transfer  coefficients  for 
Holden's  test  case. 
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TWO-DIMENSIONAL  PN$  WITH  FINITE  RATE  AND  EQUILIBRIUM 
REACTION  RATE  CHEMISTRY  ( BURROWS-KURKOV ) . 

TURBULENT  FLOW  REGIME. 

DIFFUSION  CONTROLLED  COMBUSTION. 


Test  section  intermediate 


Test  section  showing  hydrogen  injection  step,  location 
of  static  pressure  ports,  and  measurement  stations. 


Figure  22.  Schematic  of  the  Burrows-Kurkov 
Experimental  Apparatus. 
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In  this  study,  finite  rate  and 

equilibrium  reaction  rate  models  were  employed  when  computing  the 
Burrows -Kurkov  test  case.  Figures  23  and  24  depict  computed 
pressure  contours  (note  the  contour  units  are  PSF)  from  the 
finite  rate  chemistry  code  (SSCPNS)  and  the  equilibrium  chemistry 
code  (SSCPNSEQ) .  From  the  contour  plots  it  is  apparent  that  the 
results  are  quite  similar  except  that  the  pressure  rise 
associated  with  the  combustion  process  is  slightly  lower  in  the 
SSCPNSEQ  solution.  This  is  attributed  to  the  fact  that  for  this 
particular  test  case  the  degree  of  reaction  in  the  flame  region 
was  approximately  95  percent,  thus  accounting  for  the  slightly 
lowei  pressure.  The  oscillations  in  the  pressure  contours  near 
the  upper  wall  were  attributed  to  poor  grid  resolution  in  this 
region. 

The  computed  static  temperature 

contours  from  the  SSCPNS  and  SSCPNSEQ  codes  are  shown  in  Figures 
25  and  26  respectively.  As  in  the  case  of  the  static  pressure 
contours,  the  computed  temperature  contours  are  quite  similar 
except  that  the  SSCPNSEQ  solution  predicts  a  slightly  lower  peak 
static  temperature.  This  is  to  be  expected  since  the  SSCPNSEQ 
code  allows  for  the  backward  chemical  reactions.  However,  the 
maximum  discrepancy  between  the  two  solutions  was  less  than  eight 
percent. 

The  computed  mole  fraction  profiles 
from  the  SSCPNS  and  SSCPNSEQ  ^odes  are  compared  with  the 
experimental  data  of  Burrows  and  Kurkov,  with  good  agreement,  in 
Figures  27  and  28  respectively.  The  computed  total  temperature 
profile  is  shown  in  Figure  29.  Again  good  agreement  was  obtained 
with  experiment. 
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TEMPERATURE  CONTOURS 


Figure  25.  Computed  temperature  contours  from  SSCPNS 


Figure  26.  Computed  temperature  contours  from  SSCPNSEQ 


MOLE  FRACTION  PROFILES 


Figure  27.  Comparison  of  computed  and  experimentally  measured  spec 
mole  fractions. 


MOLE  FRACTION  PROFILES 

SSCPNSEQ 

(equilibrium  Chemistry) 
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Figure  28.  Comparison  of  computed  and  experimentally  measured 
species  mole  fractions. 


SUPERSONIC  COMBUSTION  SIMULATION 

8URROWS-KURKOV  TEST  CASE 
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NORMAL  DISTANCE  (CM) 


The  SSCPNS  and  SSCPNSEQ  computed 
solutions  achieved  good  agreement  with  the  data  of  Burrows  and 
Kurkov.  The  primary  difference  between  the  SSCPNS  and  SSCPNSEQ 
computer  codes  is  computational  efficiency.  The  SSCPNS  computer 
code  performed  the  above  computation  in  approximately  5  seconds 
on  a  CRAY  X-MP  computer  while  the  SSCPNSEQ  code  required 
approximately  35  seconds  to  compute  the  calculation.  The 
SSCPNSEQ  reaction  rate  model  is  non-vectorizable  and  performs 
numerous  table  look-up  operations  when  computing  the  degree  of 
reaction.  The  SSCPNSEQ  code  should  be  used  only  when  the  peak 
temperatures  predicted  by  the  SSCPNS  code  are  sufficiently  high 
to  cause  shifts  in  chemical  equilibrium. 

The  SSCPNS  computer  code  has 

demonstrated  capability  to  compute  hypersonic,  chemically 
reacting  flowfields.  The  SSCPNS  computer  code  has  demonstrated 
good  agreement  with  experimental  data  and  the  computations  of 
others  operating  at  flight  conditions. 

The  SSCPNS  formulations  allow  the 
above  computations  to  be  performed  with  little  or  no  added 
dissipation.  The  SSCPNS  code  is  ideally  suited  for  vector 
processing,  and  has  the  potential  for  optimal  parallel  processing 
acceleration.  The  SSCPNS  code  performed  the  above  computations 
accurately  and  economically,  therefore,  the  SSCPNS  code  must  be 
considered  as  a  viable  approach  to  computing  hypersonic, 
chemically  reacting  flowfields. 

This  effort  represents  a  preliminary 
study  of  the  supersonic/hypersonic  combustion  phenomena  in 
scram jet  engines. 

3.  SUMMARY 

This  document  reports  on  the  computation  of  Ramjet  Internal 
Flowfields.  It  is  divided  into  two  sections,  i.e.  subsonic  and 
supersonic  combustion.  To  numerically  simulate  subsonic 
combustion  the  time-dependent  Navier-Stokes  equations  are  used 
with  additional  species  equations  incorporated  to  model 
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hydrocarbon-air  combustion i  To  compute  supersonic  combustion  the 
parabolized  Navier-Stokes  equations  (PNS)  are  utilized  with 
species  equations  to  simulate  hydrogen-air  combustion.  Model 
cases  are  computed  and  compared  with  limited  experimental  data. 
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APPENDIX 

EQUILIBRIUM  COMPUTATIONAL  PROCEDURES 

In  the  aerospace  application  of  combustion  in  advanced 
propulsion  devices,  a  gaseous  fuel  jet  is  burned  as  it  issues 
into  air  through  a  narrow  injection  nozzle  at  high  velocity. 
Combustion  of  a  fuel  jet  is  known  to  be  a  diffusion  limited 
process.  The  fuel  in  the  jet  and  the  inlet  air  are  transported 
toward  one  another  by  convective  diffusion.  At  all  locations 
where  the  fuel  and  oxygen  are  in  stoichiometric  proportions, 
combustion  takes  place  rapidly.  In  an  ideal  diffusion  flame,  the 
reaction  zone  is  so  thin  that  it  could  be  considered  as  a  surface 
of  zero  thickness  impermeable  to  oxygen  as  well  as  fuel.  In  a 
hydrocarbon  fuel  flame  the  reaction  zone  is  visibly  blue  (flame 
temperature  about  2,000*K).  As  the  carbon  formed  in  the 
combustion  process  flows  out  of  the  reaction  zone,  it  is  usually 
incandescent  (yellow) .  Far  out  of  the  flame,  the  soot  cools  down 
to  lose  incandescence.  This  feature  is  borne  out  in  the  movies 
taken  of  the  Aero  Propulsion  Laboratory  experimental  test  results 
of  the  dump  combustor. 

To  mathematically  model  the  combustion  process  of  a 
diffusion  flame,  a  chemically  reacting  mixture  of  gases  must  be 
addressed.  When  a  mixture  of  reactants  is  placed  in  an  enclosure 
and  allowed  to  react,  surprisingly  not  all  reactants  get  consumed 
to  yield  products.  This  state  in  which  reactions  refrain  from 
proceeding  further  is  called  chemical  equilibrium.  The 
composition  depends  upon  temperature  and  pressure.  The  chemical 
equilibrium  state  is  dynamic  in  nature  rather  than  static,  i.e., 
the  reactions  never  cease  to  occur  even  at  equilibrium.  At 
equilibrium  the  rates  of  the  forward  and  backward  reactions  are 
equal . 

The  computation  of  the  equilibrium  state  involves  a  series 
of  algebraic  equations  representing  the  conservation  of  elements 
and  relationships  for  the  equilibrium  constants.  An  iteration 
scheme  is  used  to  solve  for  the  single  positive  real  root 
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produced  by  these  algebraic  relationships.  Such  a  procedure 
would  be  impractical  for  time  dependent  calculations  on  an 
unsteady  combustion  simulation.  Therefore,  an  alternate  method 
is  required,  convenient  for  approximate  simple  reactions.  The 
method  adopted  here  is  called  the  “extent  of  reaction"  method. 
It  is  described  herein  for  hydrogen/air  combustion. 

Extent  of  Reaction  Method 

Consider  stoichiometric  hydrogen-air  reaction  at  any 
intermediate  stage  of  completion,  expressed  by  the  fraction  X. 

[H2]  +  ^[02+3.76N2]  -  (1-X)  [H2+|o2]  +  X[H20]  +  ^[3.76^] 

This  equation  implicitly  satisfies  mass  conservation.  The 
information  for  this  reaction  is  tabulated  in  the  following 
table : 


i 

n . 

l 

Xi 

Y. 

l 

H2 

1-X 

U1=X1 

6.76-X 

. 028 (1-X) 

°2 

.5  (1-X) 

-OrXj. 

6.76-X 

.  226 ( 1-X ) 

H2° 

X 

2X 

6.76-X 

.  255X 

N2 

1.88 

3.76 

6.76-X 

.745 

The  single  parameter  X  may  be  determined  from  one  equilibrium 

constant  Ku  which  is  a  strong  function  of  temperature  and 
H2° 


weakly  dependent  upon  pressure. 

nH  0  ^ 

*  Vm  ■  vX 


X  -/6.760-X 
(1-X)1'5 


Hence,  X  may  be  determined  directly  from  the  thermodynamic  state 
variables  of  pressure  and  temperature. 

The  following  table  displays  the  mass  fraction  equilibrium 
values  for  the  complete  reaction  of  hydrogen-air  versus  the 
single  equation  model  (extent  of  reaction  method)  for  various 
temperatures  at  a  pressure  of  one  atmosphere. 
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TABLE  OF  MASS  FRACTIONS 


Y1 

\T 
i  \ 

2000  * K 

2  500  *  K 

3000  °  K 

3500  °  K 

p=l  atoms 

H2 

0 

0.002 

0.006 

0.009 

°2 

0.001 

0.009 

0.027 

0.038 

Complete 

Reaction 

H2° 

0.252 

0.232 

0.167 

0.063 

Model 

H 

0 

0 

0.002 

0.009 

0 

0 

0 

0.011 

0.062 

OH 

0 

0.008 

0.033 

0.056 

N2 

0.745 

0.743 

0.737 

0.729 

NO 

0 

0.005 

0.017 

0.034 

X 

0.9922 

0.9431 

0.8010 

0.5747 

H, 

0 

0.002 

0.006 

0.012 

Single 

Equation 

si  <N 

o 

0.002 

0.013 

0.045 

0.096 

Model 

H2° 

0.253 

0.240 

0.204 

0.147 

H 

0 

0 

0 

0 

0 

0 

0 

0 

0 

OH 

0 

0 

0 

0 

N2 

■  0  /45 

0.745 

0.745 

0.745 

NO 

0 

0 

0 

0 

X1 

0.992 

0.955 

0.861 

0.752 

0.0354 

0.170 

0.440 

Two 

Equation 

H2 

|  0 
° 

0.002 

0.006 

0.012 

Model 

°2 

J  0.002 

0.010 

0.031 

0.056 

h2° 

0.2  53 

C  .  235 

0 . 182 

0.107 

H 

0 

0 

0 

0 

0 

!  0 

0 

0 

0 

OH 

1  o 

0.008 

0.035 

0.080 

N2 

0.745 

0.745 

0.745 

0.745 

NO 

1  0 

0 

0 

0 

_ 
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From  the  Table  of  Mass  Fractions,  fair  agreement  is 
observed  (±4%)  up  to  about  3000°K.  Improvement  in  the  model  may 
be  achieved  by  including  additional  species.  However,  each 
species  requires  one  additional  parameter,  X It  may  be 
observed  from  the  Table  of  Mass  Fractions  that  the  next  species 
to  include  is  OH.  The  following  is  the  modification  of  the 
extent  of  reaction  method  for  a  two  equation  model. 

To  the  original  reaction  an  additional  reaction  equation  is 

added . 


Reaction  1 

[H2]  +  ^[02  +  3.76N2]  -  (1-X1)  [H^]  +  X^H.,0]  +  J[3.76N?] 

Reaction  2 

[H20]  -  *2[2H2+OH]  +  {1~X2)  [H2°] 


Combining  the  two  reactions  produces  the  following  equation: 
[H2]  +  2[°2]  -  (-1'Xl+2XlX2)  [H2]  +  2(1'A1)  [°2] 


+  X  ^ ( 1-X  2 )  [H20]  +  X1X2fOH] 


Two  equilibrium  constants  are  needed  to  evaluate  X^  and  X2 
Jp  K»  o(T) 

‘‘H2  V1,02  C1_; 

J2  XxX2 


h2° 


nH2°  ^  X1(l-X2) [6.76-X1(l-X2) ]1/2 

______ 


[1-X1+^X1X2]  [1-X1]1/2 


n 


OH 


0H(  }  7nH^  Jn0^  l1~Xi+2XlX2]1/2  C x 3 1/2 


The  information  for  this  reaction  is  tabulated  in  the  following 
table : 


V 


J 
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i 

n . 

l 

Y. 

l 

H2 

0.0283  (1-X1+  2^i^2) 

°2 

2  (1‘V 

0.226  (1-^) 

H2° 

h.  (1'A2) 

0.255  A  (1-A  ) 

N2 

1.88 

0.745 

OH 

V2 

0.241  (A1A2) 

A  comparison  of  the  two  equation  model  with  the  one 
equation  model  and  the  complete  reaction  model  is  shown  in  the 
Table  of  Mass  Fractions.  The  two  equation  model  cuts  the  mass 
fraction  error  in  half  to  about  +2%  at  3000 “K.  Property  errors 
of  this  magnitude  are  probably  satisfactory  for  numerical 
simulation,  however,  a  more  serious  difficulty  exists.  The  mass 
fractions  are  used  to  establish  the  specific  heats  (to  relate 
temperature  and  enthalpy)  and  to  ascertain  the  total  heat  of 
formation  for  the  reaction.  A  Table  of  Enthalpy  permits 
comparison  of  the  approximate  models  with  the  complete  reaction. 

The  following  is  a  summary  of  enthalpy  errors. 

ENTHALPY  ERRORS 


TfK) 

One  Equation 
Model 

Two  Equation 
Model 

2000 

0.2% 

0.2% 

2500 

0.3% 

0.2% 

3000 

1.9% 

0.9% 

3500 

6-5% 

4 . 5s- 
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TABLE  OF  ENTHALPY 


P=1  atmos. 

T  (  ‘  K) 

Yihi 

Molecule 

2000 

2500 

3000 

3500 

H2 

0 

30 

130 

230 

°2 

10 

170 

630 

1070 

H2° 

4380 

5490 

5040 

2330 

H 

0 

0 

30 

140 

0 

0 

0 

150 

990 

OH 

0 

140 

700 

1450 

N2 

10000 

13200 

16330 

19400 

NO 

0 

90 

390 

920 

h (cal) 

14390 

19120 

23400 

26530 

H2 

0 

30 

130 

300 

°2 

30 

240 

1060 

2710 

H2° 

4390 

5680 

6160 

5430 

H 

j  0 

0 

0 

0 

0 

0 

0 

0 

0 

OH 

0 

0 

0 

0 

N2 

10000 

13230 

16500 

19820 

NO 

0 

0 

0 

0 

h (cal) 

14420 

19180 

23850 

28260 

H2 

0 

30 

130 

300 

°2 

30 

190 

730 

1580 

H2° 

4390 

5560 

5500 

3950 

H 

0 

0 

0 

0 

0 

0 

0 

0 

0 

OH 

0 

140 

750 

2070 

N2 

10000 

13230 

16500 

19820 

NO 

0 

0 

0 

0 

h (cal) 

14420 

19150 

23600 

27720 
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Again,  it  is  apparent  that  the  heat  capacity  of  the  mixture  can 
be  satisfactorily  produced  up  to  3000 *K  with  either  model.  This 
is  not  true  for  the  heat  of  formation,  however.  A  comparison 
with  the  Table  of  Heat  of  Formation  shows  the  growth  of  large 
errors  with  both  models  as  shown  below: 

HEAT  OF  FORMATION  ERRORS 


Tr  *ki 

One  Equation 
Model 

Two  Equation 
Model 

2000 

0.7% 

0.7% 

2500 

6.0% 

3.0% 

3000 

85% 

60% 

3500 

204% 

165% 

Clearly,  an  improvement  is  required  in  order  to  obtain 
correct  values  for  the  heats  of  formation.  Since  the  one 
equation  model  produces  satisfactory  mass  fractions  and  specific 
heats,  it  is  proposed  that  this  one  equation  model  be  retained 
with  modified  heats  of  formation.  The  heat  of  formation  for  H20 
is  -60.5  kcal/mole  for  the  range  of  interest.  The  formation  of 
trace  species  greatly  reduces  the  heat  of  combustion  and,  in 
fact,  produces  positive  values  at  3500°K  due  to  dissociation 
effects.  The  following  is  a  table  of  the  heats  of  combustion  for 
stoichiometric  hydrogen-air  for  the  complete  reaction  model  at 
three  different  pressures. 


i 
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TABLE  <5F  HEAT  OF  FORMATION 


- - 1 

P=1  atmos. 

T  ( 0  K) 

Xihi 

Species 

2000 

2500 

3000 

3500 

H2 

0 

0 

0 

0 

Complete 

°2 

0 

0 

0 

0 

Reaction 

H2° 

-20630 

-18770 

-12590 

-  4070 

Model 

H 

0 

+  220 

+  2250 

+  9500 

0 

0 

+  60 

+  980 

+  4500 

OH 

+  8 

+  90 

+  360 

+  500 

N2 

0 

0 

0 

0 

NO 

4-  11 

+  90 

+  280 

+  470 

Ah, (cal) 

-20610 

-18300 

-  8720 

+10900 

t 

H2 

0 

0 

0 

0 

Single 

°2 

0 

0 

0 

0 

Equation 

H2° 

-20750 

-19400 

-16160 

-11350 

Model 

H 

0 

0 

0 

0 

0 

0 

0 

0 

0 

OH 

0 

0 

0 

0 

N2 

0 

0 

0 

0 

NO 

0 

0 

0 

0 

Ah, (cal) 

-20750 

-19400 

-16160 

-11350 

£ 

H2 

0 

0 

0 

0 

°2 

0 

0 

0 

0 

h2° 

-20750 

-18950 

-14280 

-  7950 

H 

0 

0 

0 

0 

O 

0 

0 

0 

0 

OH 

0 

+  100 

+  300 

+  850 

N2 

0 

0 

0 

0 

NO 

0 

0 

0 

0 

Ah, (cal) 

-20750 

-18850 

-13980 

-  7100 
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HEATS  OF  COMBUSTION 


t 


i 


Ahf 

kcal/mole  of 

mixture 

'"^\T(*K) 

P  ( atoms}'\^^ 

2000 

2500 

3000 

3500 

in 

• 

o 

1 

-17.58 

-  5.22 

1.0  1 

-18.31 

-  8.72 

+  10.9 

2 . 0 

-18.94 

-11.50 

+  5.06 

V 

in 

• 

o 

0.302 

0.178 

1.0 

0.343 

0.311 

0.208 

0.067 

2.0 

0.318 

0.234 

0.0995 

Ah  £ 

kcal/mole  H2 

0 

0.5 

-60.5 

-58.21 

-29.32 

1.0 

-60.5 

-58.87 

-41.92 

+  163 

2.0 

-60.5 

-59.56 

-49.14 

+  50.8 

By  using  a  reduced  value  for  the  heat  of  formation  as  shown  in 
the  above  table,  all  thermodynamics  of  the  combustion  flow  may  be 
simply  modeled. 
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